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Abstract 

Based on a dipolar-elastic model for oxygen vacancies on rutile (110), we evaluated analytically 
the overall energy of a periodic array of two vacancies and extracted the interaction parameters 
from total-energy density functional theory (DFT) calculations. Our calculations show that the 
dipole model holds for next-nearest neighbor vacancies and beyond. The elastic-dipolar interaction 
vanishes for adjacent vacancies, but they still experience an electrostatic repulsion. The proposed 
interaction model predicts a vacancy separation distribution that agrees well with that determined 
in our ultra-high vacuum scanning tunneling microscopy experiments, and provides a perspective 
for understanding earlier DFT reports. 
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Titanium dioxide -widely used in heterogeneous catalysis [1], photocatalysis [2[, solar 



cells 



3|, or gas sensors 



of metal oxide surfaces 



4 1, has become the prototype material for studying the reactivity 
5[. Defects such as oxygen vacancies are always present on rutile 
surfaces fl and, depending on their coverage and spatial distribution, can strongly influence 
the reactivity of the surface [7]. The interactions between vacancies determine their spatial 
distribution on the surface. Highly reactive vacancy clusters or pairs have not been expected 
to form because of vacancy repulsions but recent experiments (9] do show the possibility 
of spontaneously formed oxygen vacancy pairs (OVPs), i.e., of two adjacent vacancies in 
the same bridge-oxygen row. Regarding the stability of OVPs, early density functional 
theory (DFT) calculations came to contradictory conclusions. The OVPs were reported to 



have the highest 10] and the lowest ll| energy of all configurations of two vacancies per 
computational cell. A newer study jg] finds virtually the same energies for the OVP and 
the next-nearest neighbor (NNN) configurations, while another recent study |q| reports the 
NNN structure to have a much higher energy than the OVP. To date, several issues have 
prevented the complete, fundamental understanding of vacancy interactions, including their 
reliable quantitative determination; the more important issues are the difficulty of decoupling 
the interactions while using computational slabs of manageable size, and the sensitivity of 
various structural properties to the number of layers in the supercells 
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Here we show that the interaction of same-row vacancies on rutile (110) is dipolar-elastic 
in nature, with a long-range, inverse-cube dependence on their separation. This dipolar- 
elastic model holds when the vacancies are not adjacent, which we have found from DFT 
calculations at the level of the generalized-gradient approximation (GGA). Our approach 
has two key features that allow us to reliably determine the formation energies and the in- 
teraction parameters from total-energy GGA calculations: first, the interactions have been 
isolated to one bridge-oxygen row by using large supercells, and second, we have developed a 
closed-form expression for the overall interaction (per computational cell) associated with a 
periodic array of two vacancies. When vacancies are adjacent, they still repel, but this repul- 
sion is much weaker than the dipolar-elastic one at the same distance. We have determined 
the distribution of vacancy separations (along bridge-oxygen rows) by scanning tunneling 
microscopy (STM), and have found that this distribution agrees well with that predicted 
from the calculated interactions. This validates our physical model for vacancy interactions, 



which we use to analyze our DFT data as well as data from other works 



BE 




FIG. 1. (Color online) Reduced 10 x 2 rutile surface slab used in the DFT calculations. The 
interaction between the two bridge-oxygen vacancies is determined for different values of their 
separation d, c < d < 5c. 

For the DFT simulations, we constructed 10 x 2, 600- atom stoichiometric supercells with 
dimensions L x = 2ay/2 and L y = 10c (with a = 4.669 A, c = 2.970 A [lsj), and a thickness 
of five O-Ti-0 trilayers. The vacancies were created by removing two same-row oxygen 
atoms spaced at d (c < d < 5c) [Fig. [T] . The DFT relaxations were carried out in the GG A 
ramework using the PBE exchange-correlation functional [jjj], projector-augmented wave 



151 ] pseudopotentials 



16[, and an on-site Hubbard term U for the Ti 3d states [17|. Charge 



neutrality (Q = Oe) was maintained for the stoichiometric slabs, but for the reduced slabs 
we also considered the positively-charged case (Q = 4e, corresponding to the removal of two 
2 ~ ions). We have not searched for the localized electron configurations that optimize the 
total energy 18|, but simply relaxed the structures from the bulk truncated positions and 
analyzed their final electronic structures. For our non-zero Hubbard term values, we have 
found that localization occurs on subsurface Ti atoms for all spacings d > c. 

The difference AE between the total energy of the reduced slab (E r ) and the energy of 
the stoichiometric one of same area and thickness (E s ) can be written as 



AE = E r -E s = 2(f-fi )+w, 



(1) 



where / denotes the formation energy of a single vacancy (on an otherwise perfect and wide 
surface), fio is the oxygen chemical potential, and w contains all interactions. Since we have 
collected all interactions into a single term w, which depends on the supercell dimensions and 
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FIG. 2. (Color online) (a,b) The difference between the energy of a reduced (10 x 2) slab [(a) 
neutral, (b) positively charged] and that of a stoichiometric one. (c) Analytical dependence u{rj) 
(Eq. ([5]), curve) and numerical calculations of u for neutral (x symbols) and charged slabs (+ 
symbols). 



on the spacing between the vacancies, the formation energy / in Eq. (JTj) depends neither 
on the spacing between vacancies nor on their coverage. The variation of AE with the 
separation d at constant L y (L y = 10c) is plotted in Figs. I2^a,b) for neutral and positively 
charged slabs. In order to extract interaction parameters from AE vs. d data, we have to 
understand the overall interaction term w. 

Therefore, we first focus on finding an analytic expression for w, and start by neglecting 
the cross-row interactions; this is reasonable given the large supercell dimension along [110], 
L x = 2ay/2 = 13.205A. In this approximation, w [Eq. (pQ)] depends only on the vacancy 
separation d and on the dimension L y along a bridge-oxygen row. In the framework of 
elasticity theory, point defects on surfaces interact as elastic multipoles whose long-range 
interactions are inversely proportional to certain powers of their separation [3]. In what 
follows, we describe the interaction v between two isolated vacancies by the long-range 
dipolar-elastic repulsion 



v(d) 



V\ if d 
§ if d 



ic, i = 2, 3, 4, 



(2) 



where d is the distance between the two vacancies on an otherwise perfect surface, G is the 
strength of the dipolar repulsion, and v± is a short-range interaction present only for adjacent 
vacancies. When using periodic boundary conditions, the two vacancies are not isolated, 
since they interact with their periodic images as well. Using ()2]) for d > c and collecting the 
contributions from all periodic images along the same row, the total interaction energy per 



supercell can be written as a function of L y and rj = d/L y via 



G 

w{L y ,r])^ —u(t]), (3) 



with the function u(rf) given by 



^ 3 + ^ U 3 + (^ + ^) 3 + (k- 



= (1/v 3 ) + 2C(3) - (^ 2) (v) + ^ (2) ("^))/2 (4) 
where ( is the Riemann zeta function (C(3) — 1.202) and zM 2 ) is the polygamma function of 



second order. Polygamma function identities 



20| reduce Eq. (0} to 



u(rj) = 2C(3) - vr 3 cot (ttt?) csc 2 (tt//) - V' (2) (77) • (5) 

Eq. (jSJ) is a general description of the interactions of two identical, elastically-repelling defects 
in the same row and their periodic images along that row; as such, it does not depend on 
the (common) type of the defects (e.g., both vacancies or both adatoms), on their formation 
energy, or on their interaction strength. 

TABLE I. Formation energies / (= / + fJ-o), repulsion strengths G, and short-range interactions 
v\ for different values of the slab charge Q and Hubbard parameter U. The standard deviations 
for / and / are the same. 



Q(e), U(eV) 


/(eV) 


/(eV) 


G(eVA 3 ) 


ui(eV) 


0, 0.0 


7.170 ±0.006 


2.244 


169.8 ±4.1 


0.677 ±0.013 


0, 3.0 


7.611 ±0.013 


2.684 


152.0 ± 10.0 


0.702 ± 0.027 


4, 0.0 


16.966 ± 0.004 


10.207 


157.7 ±3.3 


0.575 ± 0.009 


4, 3.0 


16.964 ± 0.006 


10.205 


165.1 ±4.6 


0.808 ± 0.013 



Using Eqs. ([T]), ((3]), and ([5]), we fit the data in Figs. [2] (a,b) for d > 2c to obtain the 
relative formation energies / = / — fio and the interaction strengths G for different Q (slab 
charge) and U (Hubbard parameter). The fio values [see Eq. (OQ)] that we have used were 
fio — —4.926 eV (half the energy of an O2 molecule) for the neutral system, and fio — —6.759 
eV (the energy of an isolated O 2- ion in the supercell) for Q = 4e. The calculated formation 
energies and interaction strengths are listed in Table [J for neutral and charged slabs at 
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FIG. 3. (Color online) Displacement fields in a plane containing the vacancies separated by (a) d = 
5c and (b) d = c; for clarity, only three trilayers are shown. The green arrows show schematically the 
horizontal force dipoles associated with each vacancy (a), and illustrate the monopole cancelation 
responsible for the vanishing elastic repulsion at d = c (b). The displacements (magnified here 
10-fold for clarity) are calculated for the (Q = Oe, U=0eV) system with respect to the relaxed 
stoichiometric slab. The largest displacement magnitude is 0.41 A in (a) and 0.44A in (b). 



U = 0.0 eV and U = 3.0 eV. Using / and G values from Table HI we check our model by 
numerically calculating u from Eq. ([[]) [i.e., u = (AE — 2f)Ly/G for rj > 2c/ L y ] and plotting 
it along with the analytic result Eq. (JSJ). As seen in Fig. |2^c), the agreement between the 
numerical u values and the general formula Eq. (jSJ) is very good, which validates a posteriori 
our assumption that the interactions are reasonably well-confined to the bridge-oxygen row 
as long as the neighboring rows are defect-free. 

Although multipole interactions between atomic-level defects (most often adatoms) on 



crystal surfaces have been studied (2l|, so far the particular dipolar-elastic model proposed 
here has not been reported for oxygen vacancies. Figure [31(a) shows the atomic displacement 
fields and, schematically, the horizontal force dipoles (F + , F_) associated with each vacancy 
for d = 5c. The atoms located between vacancies experience opposite pulls resulting in an 
increase of energy, i.e., the elastic repulsion. When the vacancies are brought close to form 
an OVP, there are no more 5-fold coordinated Ti atoms (5-f Ti) between them, which leads 
to the cancelation of two force monopoles as shown in Fig. [3](b). It may be worth noting 
that monopole cancelation has also been reported to be the origin of a short-range attraction 
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that leads to step bunching on certain surfaces [22j. Despite this monopole cancelation, the 
interaction between vacancies at d — c is not zero but a quantity v\ [Eq. (J2])] that can be 
found from a straightforward modification of Eq. (131) , 

c G G 
w{L y ,r) 1 = —) = —u{r)i)- — + v 1 . (6) 

±Jy lJ y C 

Using Eqs. ([I]), §6§ and the / and G values already calculated, we have found that 
the interaction v\ is small but positive in all cases [Table HJ last column]. This short- 
range repulsion is about one order of magnitude smaller than what the dipolar-elastic model 
would predict for vacancies at d = c [G/c 3 ~ 6.48 eV], and is largely due to the electrostatic 
interactions of the exposed 4-f Ti with the nearby 5-f Ti atoms. 

The directly observable manifestation of vacancy interactions is their spatial distribution, 
which we have analyzed from thermal vacancy populations. The samples were produced by 
Ar + sputtering cycles, followed by annealing the (110) rutile surface at a temperature T = 
950 K, then rapidly cooling down to 77 K to freeze in the vacancy distribution. The vacancies 
were imaged using an ultrahigh vacuum (UHV, pressure below 10~ n Torr) cryogenic STM, 
in constant current mode with positive sample bias {23 1. Vacancy separations were analyzed 
from portions of bridge-oxygen rows that had a vacancy concentration of n ps 15%. If 
there were no interactions, then a fixed vacancy coverage n would lead to an exponential 
decay of the probability to find vacancy-to- vacancy (V-V) segments of length d, p n onint(d) oc 
exp (— nd/c) {2^]. Our data shows that p(d) exhibits a maximum, and not the monotonic 
decay corresponding to the non- interacting system [green curve in Fig. [4]; this is a direct 
consequence of the repulsive interactions between vacancies. 

In a canonical ensemble system of V-V segments, the interactions between the ends of 
segments give the single-particle energy levels v(d) (d = c, 2c, 3c, ...) and thus a canonical 
distribution p(d) = {l/Z) exp (—nd/c) exp (—v(d) /ksT), in which Z is a normalization factor 
and ks is Boltzmann's constant. The canonical distribution based on the interaction model 
Eq. (J2J) with the parameters in Table |T] is consistent with the experimental data [refer to 
Fig. H], and is virtually the same for all (Q, U) pairs used in this study. The competition 
between the fixed coverage constraint and the rapidly-decreasing dipolar repulsion (G/d 3 ) 
gives rise to a most-probable vacancy spacing d* which can be readily derived from the 
canonical distribution, d* = (?>Gc/nkBT) l / A . The G values (Tabled]) give d* between 6.1c 
and 6.3c, consistent with the experimental peak at 6c. The agreement between the vacancy 




FIG. 4. (Color online) Distribution of vacancy spacings for n = 15% vacancy coverage at 950 K: 
STM experiments (dots with error bars) compared with the non-interacting system (exponential 
decay, green curve) and with the canonical distribution (four nearly overlapping curves) . Note that 
the canonical distributions for different Q and U values are not fits to the experimental data, but 
are based on the interaction model Eq. ([2]) with the GGA values for G and v\ listed in Tabled 



separation statistics determined in STM and the canonical distribution with GGA-calculated 
parameters validates our interaction model Eq. (J2J). 

In the previous systematic attempts to compute the vacancy interactions on TiC>2(110) 
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11], the total energy was expressed as pairwise interactions from each vacancy to the 



id, 



llj acknowl- 



next one along the row and also included cross-row couplings. Both reports 
edged unresolved shortcomings of the pairwise model, which had manifested in significant 
differences of the total energies predicted by the model (with respect to those obtained di- 
rectly from DFT calculations) once the model was applied to supercells other than those 
used to determine the pairwise interaction parameters. Departing from these pairwise mod- 
els of Refs. IOj, |ll(, we have proposed herein that the same- row vacancy interactions are 
dipolar, thus long-ranged. As we will show below, our model Eq. ([H])-© holds very well 
when applied to different supercells (than those in Fig. [1]), different numbers of vacancies 
(one or two) per cell, and different exchange-correlation functionals. For example, we have 
used Eqs. flSl)-® and the (Q — 0, U — 0) case data in Table fl] to compute the total energy 
difference between 5x3 supercells with two vacancies in the NNN and OVP configurations. 
We have obtained 0.264 eV, in excellent agreement with our GGA simulations performed 
for 5 x 3 slabs, which place the NNN supercell energy at 0.249 eV above that of the OVP 
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supercell; for the other three cases in Table HJ the model [Eqs. ©-(JS])] yields total energy 
differences that are within 0.06 eV or less from the GGA results. Our results for 5x3 super- 
cells are in quantitative disagreement with the recently reported total energy difference of 



0.8 eV js], but based on our calculations and on the STM experiments that show similarly 
small occurrence probability for OVPs and NNNs (Fig. H]), we believe the 0.8 eV value to 
be in error. 

Interestingly, data from other reports of DFT simulations can be readily understood 



using our model. The two-vacancy results in Ref 



.a 



(c < d < 5c) can be fitted well to our 



Eqs. ©, ©, yielding a strength G = 244.9 ± 7.9 eVA 3 (not enough data is provided to 
determine / or /). For one vacancy per supercell, the quantity denoted as VFE in ljj is 
defined as our AE = E r — E s up to an additive /xo- We find that the data for p(m x 1) 
cells (m = 2,3,4,6) in table 2 of Ref. [ll] fits closely our dipolar-elastic model, which for 
single-vacancy supercells takes the simple form AE + /xo = / + G£(3)/Ly. This data yields 
/ = 3.134 ±0.012 eV and G = 198.98 ±4.06 eVA 3 ; these values are consistent with those in 
the first line of our Table [J but differ from them likely because of the different computational 
parameters used in While we have devised our physical model for interactions confined 



to the same row, the structures in Refs. 



11] do not have any intact oxygen row and 



thus allow for cross-row coupling. Even so, the p(m x 1) data fits our elastic model very 
well, as judged by the small standard deviations obtained for / and G; in the case of 



Ref. [HI] (p(m x 1) cells in table 2, m > 2), this is because at one vacancy per supercell, the 
cross-row interactions occur mostly perpendicular to the oxygen rows and thus may amount 
to a constant independent of m (the cell dimension along the row). For the two- vacancy 
results in Ref. , the agreement with our model likely occurs because the diagonal cross-row 
interactions do not vary significantly as a function of d when d > c. 

In conclusion, we have shown that the dipolar-elastic model describes well the long-range 
repulsion of same-row vacancies for all separations except d = c, where a much smaller short- 
range interaction is present. The model fits not only our DFT data, but also explains several 
other results from the literature and gives an equilibrium vacancy separation distribution 
that agrees well with that determined in our STM experiments. 
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